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Abstract. 

In this work we introduce a differential equation model with time-delay that 
describes the three-stage dynamics and the two time scales observed in HIV infection. 
Assuming that the virus has high mutation and rapid reproduction rates that stress 
the immune system throughout the successive activation of new responses to new 
undetectable strains, the delay term describes the time interval necessary to mount 
new specific immune responses. This single term increases the number of possible 
solutions and changes the phase space dynamics if compared to the model without 
time delay. We observe very slow transits near the unstable fixed point, corresponding 
to a healthy state, and long time decay to the stable fixed point that corresponds to the 
infected state. In contrast to the results obtained for models using regular ODE, which 
only allow for partial descriptions of the course of the infection, our model describes 
the entire course of infection observed in infected patients: the primary infection, the 
latency period and the onset of acquired immunodeficiency syndrome (AIDS). The 
model also describes other scenarios, such as the very fast progression to the disease 
and the less common outcome in which, although the patient is exposed to HIV, he/she 
does not develop the disease. 
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1. Introduction 

After almost three decades of research attempting to understand the dynamics of HIV 
infection, many aspects of its underlying mechanisms remain unrevealed. Understanding 
the three-stages dynamics and two time scales has also challenged the scientific 
community over the past decades. Despite the intrinsic differences between individuals, 
the evolution of HIV infection follows a common pattern [1] that starts with the primary 
infection during the first weeks after contamination and is characterized by an extensive 
dissemination of the infection (large virus count), followed by a pronounced decline 
in the virus count caused by the development of the virus-specific immune response. 
Notwithstanding the decrease in the virus load, the system does not recover completely 
from the infection after the primary infection and very low concentrations of virus 
remain in the organism. This second stage, called the latency period, is asymptomatic 
and varies from patient to patient. The time-scale of the clinical latency period ranges 
from months to several years, and if untreated, is characterized by a progressive decrease 
in the number of HIV target cells, the TCD4 + cells. The third phase corresponds to 
the onset of acquired immunodeficiency syndrome (AIDS) defined as the time when the 
T cell counts reaches the order of 20 — 30% of the concentration of the TCD4 + cells 
in healthy individuals [2]. Without any treatment, the patient dies from opportunistic 
diseases. 

Among the efforts to understand the dynamics of the disease, many mathematical 
models have been proposed to describe either specific aspects of the HIV dynamics or 
the overall behavior of the evolution of the infection. The attempts of the first two 
decades are reviewed in details in [3111] and references therein, but recent contributions 
still indicate an active interest in the theme [5H?]. To date, most of the proposed models 
have been based on differential equations (ODE and PDE) ( see for exemple [3^l8]-[T2"] ) . 
although discrete models have also been considered [TBTTTd] . While most of the models 
proposed until 2001 were very successful in describing specific aspects of the dynamics or 
the primary infection or the latency period, none was able to describe the entire course of 
the infection employing the same set of parameters. The first model that described the 
three-stage dynamics and two time scale as observed in infected patients was a cellular 
automaton model put forward to describe the dissemination of infection in the target T 
cells located in the lymph nodes [15] . Recently, Stilianakis and Schenzle [XT] proposed an 
ODE approach to reproduce the entire course of HIV infection. The model describes the 
targeted infection of TCDA + cells generating new variants of HIV that compete between 
themselves with a capacity of increased reproducibility within the selected variants that 
in turn leads to the immune system's loss of control over the infection. Using seven 
non-linear ODE' s and 22 parameters, the results presented by the authors reproduce 
the entire course of the infection on the time scale of days, corresponding to the very fast 
progression of the disease. The authors claim that although the results presented do not 
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reproduce the typical course, they were able to reproduce (not revealed in the paper) 
the whole spectrum of courses observed in infected patients by varying the parameters. 
In a further discussion, they mention that due to complications in the dynamics of the 
disease, the model loses its applicability to the third phase of the infection. According 
to the results and discussions presented throughtout the paper the model seems to 
reproduce the three stages only for fast progressions of the disease (one time scale) and 
fails to reproduce the two time scales (days and years) observed in infected patients to 
describe the third phase. 

Here we seek to describe the entire course of the HIV infection using a unique 
set of parameters and a non-linear differential equation approach that includes time- 
delayed terms. Our model follows the same assumptions as the cellular automata 
model introduced by Zorzenon dos Santos and Coutinho [15]: HIV's high mutation 
and rapid replication rates and the ability of the immune system to mount a new 
regular immune response to any non-identified HIV resulting from the virus proliferation 
process. However, since a time interval (r) is necessary for the immune system to mount 
any specific response developing from the necessary signaling among the cells, the non- 
identified strain remaining during this time interval is able to disseminate the infection 
throughout the target cells. In the follow up, these newly infected cells would be able 
to produce new strains (through mutations) that may not be recognized by the immune 
system (IS) therefore, disseminating the infection during the same period of time. 

The addition of this single delay term changes the dynamics of the model in the 
phase space, generating trajectories that depend on its non-local properties. These 
properties change the transit of the trajectories near the fixed points allowing for the 
description of the entire course of the HIV with the two time scale and three-stage 
dynamics. Time -delay effects are present in different (molecular, cellular or organ) 
mechanisms of various biological systems [18] and recently, other models with time 
delay terms have been introduced to describe HIV infection. However, in these models, 
the time delay either represents a finite incubation time before the release of new viruses 
by infected TCDA + cells [19,20] or describes the delayed effects of the anti-retroviral 
drugs on the plasma viremia that leads to the viremia decay [3 [11] [ZD 122] • Although 
these ODE models have time delay terms, they do not reproduce the two time scales 
and three-stage dynamics of the HIV infection as in the model we introduce. 

2. The model 

In order to facilitate the understanding of the model introduced herein, let us briefly 
summarize the main assumptions of the CA model |15j . The CA model describes the 
local interaction among target cells (T cells and macrophages) and the virus in the 
lymph nodes, taking into account the HIV's high mutation and rapid reproduction 
rates, as well as the fact that the immune system responds regularly to HIV as to any 
other virus. The target cells can be found in four states: healthy (h), infected (A and 
B) and dead (d) (or empty sites). Healthy cells are the target cells that may become 



The dynamics of the HIV infection: a time-delay differential equation approach 4 



infected. Infected-A corresponds to the stage in which the infected cell disseminates the 
infection during r time steps (the period of time necessary to mount the specific immune 
response) and infected-B corresponds to the last state of the infection in which the cells 
are recognized and killed in the next time step. The dynamics of the interaction among 
automata is defined by four rules: (1) h cells become infected if at least one nearest 
neighbor is infected- A. Since it is considered that each new infected- A that enters 
the system carries a new strain of virus, (2) asserts that infected -A cells spread the 
infection into its vicinity during r time steps, and afterwards the infected-A cell turns 
to a less infective-B stage. Rule (3) describes the turnover of B cells into d cells (or 
empty sites) in the next time step. The regular blood flow in the system allows for the 
natural replenishment of the target cells in the lymph nodes, therefore (4) replaces d 
cells by h with probability p rep i or by infected-A cells with probability p rep i * Pi n f ec , 
since we also consider that the empty sites may be occupied by infected cells coming 
from other lymphatic compartments. In the CA model, the fast dynamics of the primary 
infection corresponds to a rapid increase in the number of infected cells through the cycle 
h— >A— »B— >d— >h. The slow dynamics of the latency period is related to the formation 
of the spatial structures of infected cells. As shown in [15], these growing structures 
slowly compromise more and more cells, segregating and trapping the healthy ones and 
leading to a reduction in TCD4 + cell counts. These structures are associated with 
syncytia (aggregation of infected cells) formation, observed in HIV cultures and in the 
lymph nodes of HIV patients. 

Our model describes the time evolution of TCD4 + cell density in healthy, infected 
and dead states using a set of ordinary differential equations. As in the CA model, 
we have differentiated the stages A and B in the population of infected cells, but in 
turn there are also two differentiated two stages in the population of healthy cells : h\ 
corresponding to the healthy cells present in the tissue at the beginning of the HIV 
infection, and h 2 corresponding to the population of new cells that enter the system 
when the infection has already set in. If instead we consider only one type of healthy 
cells we will obtain the same qualitative results, but shorter latency periods, as discussed 
bellow. The sum of all the variables (populations) remains constant over time and the 
time delay terms would be included in the equations for infected-A and -B cells to 
describe the time necessary for the IS to convert any new A cell into a type B cell. 

Therefore, the model is described by the following set of differential equations: 

h x = - k 5 h (t)A p - k 6 hB n , 
h 2 = k 3 d - hh 2 A q - k 6 h 2 B n , 

A = — k\A (t — t) + hd + hihAP + h 2 A q ) + k 6 (h + h 2 )W\ 
B = k x A (t - r) - k 2 B (t) , 

d = -k 3 d(t)-hd(t) + k 2 B(t). (1) 

where p and q represent the number of A neighbors required to convert healthy (hi 
or h 2 ) cells into newly infected-A cells; n is the number of B neighbors necessary to 
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infect healthy cells (hi or h 2 ). As in [J5], we consider p = 1, n = 4 and r = 4. The hi 
population has no correspondent in the CA model but is expected to be more resistant 
to HIV, requiring a larger number q (when compared to p) of A neighbors to become 
infected. This assumption is based on the fact that once the virus is detected and the 
specific immune response is developed the system would respond faster to newly infected 
cells, therefore more infected cells would be necessary to infected a healthy one, once 
the infection has set in. The rate constants describe the different transitions between 
states. Two of them can be directly related to CA parameters: k 3 being the rate at 
which new A-cells are produced is related to p rep i mentioned above, and /s 4 corresponds 
to the rate at which newly infected cells enter the system( p rep i *Pi n f ec )- The remaining 
rate constants, fci, k 2 , k§ and k§ describe, respectively, the transitions from states A to 
the B, B to d, and hi and h 2 to A. 

As in the CA model [15] , values were chosen for some rate constants and parameters 
based on certain biological data : (i) the replenishment ability of the system k 3 0(1) ; 
(ii) a very low but finite probability (1 in 10000 cells of the peripheral blood harbor the 
viral DNA) that some of the newly infected cells entering the system would come from 
other compartments hi « 1; and (iii) since the responses to different virus can vary 
from a few days to 8 weeks, as in the CA model, we adopted r = 4. 

The values of ki are not freely chosen, but obey certain constraints that result from 
considering : i) the coordinates of the fixed points (FP = (hi, h 2 , A, B, d) representing 
densities that should remain in the [0, 1] interval; (ii) at the infected state the densities 
would respect the inequalities: infected — A > infected — B > and d < infected — B 
; and (iii)p rep ; > p in f otherwise would replace the empty sites only by infected-A cells. 
Therefore, k 2 > ki, k 5 > ki, k 3 + A; 4 > k 2 , k 3 > fc 4 . 



3. Results 



The system of equations (CQ) admits three FP solutions. Two, FP = (1,0,0,0,0) and 
FP = (0,1, 0, 0, 0), represent the healthy states and a third that represents the infected 
state FP\ = (0,h 2 ,A,B,d). To find its coordinates, it is necessary to solve the following 
equation for B : 



RB + SB 



T 



k 2 

h 



9-1 



B q + U 



B*- 1 + V 



0. 



(2) 



Equation ([2]) is obtained by making the left-hand sides from the equations of system 
(GQ) equal zero, as well as considering the condition of normalization. Once equation 
(J2J) has been obtained, it is then possible to describe the behavior of the infected B- 
cells, in which coefficients R, S, T, U and V are functions of the parameters of the 
model. Depending on the value of q ( integer or not), (J2J) becomes either a polynomial 
or a transcendental equation. However, for integer q e [0,4], the resulting analytical 
expression for B is very complex and has no practical application to the analysis of FP\. 
Therefore, (j2J)is always solved by numerical methods. 
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Figure 1. Location of some eigenvalues in the (a, /3) plane, corresponding to the 
intersections of the nullclines u(a, f3) = (solid) and v(a,f3) = (dotted). The main 
features are: u = when fj = 0; there are two real eigenvalues for a — (Ao), and 
a < (A4);there is a pair of complex eigenvalues (^2,3) close to the origin; and an 
absence of eigenvalues with a > 0. 

FP and FP are both unstable. FP is reached only for initial condition with 
hi=l, while few heteroclinic orbits converge to FPq along the direction of its attractive 
manifold. To analyze the stability of the FP%, we follow the standard procedure using 
the equation for the eigenvalues of the Jacobean matrix J, M(X)=det(J — AJ)=0. For 
r >0, M(A)=0 becomes a transcendental equation since it contains terms that depend on 
e~ Ar . Irrespective of the value of r, we always obtain an identically vanishing eigenvalue 
(Ao), related to the conservation of the number of cells. 

When r=0 there is no delay in mounting the specific immune response 
and its spectral equation is polynomial. For example, if we choose k\=0.l6, 
k 2 =0.2, /c 3 =0.68, k A =2.5 x 10~ 5 , k 5 =0.6, k 6 =0.1, n=4, p=q=l, we obtain FPi=(0, 
0.266,0.361,0.289,0.085), while the spectrum is given by (A , A l5 A 2 , 3 , A 4 ) = (0,-0.705,- 
0.196±0.215i,-0.217). Since hi = 0, A4 describes only the decay of hi to 0; Ai indicates 
a fast decay to a two-dimensional space in which the slow decay towards FPi , actually 
described by A2,3, takes place. In other words, as observed in the CA model [23J the 
absence of time delay in mounting the specific response would favor the spread of the 
infection and shorten the latency period. If instead of unity we adopt q = 1.13, the 
decaying dynamics slow down as Re{\2 t z) — —0.169, i. e., the effect of shortening the 
latency period can be reduced if we increase the resistance of the new incoming healthy 
cells. 

For r > 0, the complex eigenvalues are obtained employing the Newton- Raphson 
(NR) procedure, using the notation A = a + i(3 where a = Re(X) and = Im(X), 
M(A) = u(a, {3)+iv(a, (3). The eigenvalues correspond to the points where the nullclines 
of u(a,/3) = and v(a,(3) = intersect in the (a, (3) plane, as illustrated in Figure 1. 
Graphs are very helpful for finding roots of M(A) = 0, since they indicate the initial 
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values for the NR procedure. The overall picture remains essentially the same if we 
enlarge the region of the (a, /3) plane in which the nullclines are drawn. Note that the 
nullclines cross each other only in the negative a region and FP\ remains stable. The 
several nullcline branches, when a < 0, indicates that both u and v oscillate and that 
there are an infinite number of solutions to M(A) = 0. The influence of r on the decay 
rate influences only the eigenvalues with the smallest real part, i.e., A2,3- The stability 
analysis shows that, when r is increased from to 4, the decay time to FPi is amplified 
by factor 7 for q = 1.13, and factor 5.4 for q = 1.0. All results of the stability analysis 
for r = are checked by measuring the slope of the amplitude of the decaying solution 
to FP\ that coincides with Re(\2^)- 

If we vary the parameters respecting the conditions discussed at the beginning of 
this section, we note that the coordinates of FPI change in a continuous manner and 
FP\ remains a stable fixed point. Furthermore, by fixing the values of the rate constants, 
we find that the attracting properties of FP X weaken as q and r increase. These two 
parameters of the model contribute to reducing the velocity at which the trajectories 
decay. In other words according to this model the duration of the latency period is 
determined by the slow dynamics on the trajectories close to the stable fixed point, 
which in turn is controlled by parameters q and r. By increasing the value of q we require 
a larger number of infected- A cells in order to infect healthy cells /12, indicating that an 
increase in the resistance of the incoming target cells also increases the latency period. 
On the other hand, greater the r, the greater the latency period or the individual' s 
survival lifetime. Nevertheless, these features are not sufficient in themselves to account 
by themselves for the observed increase in the time-scale associated with the latency 
phase. It also depends on the presence of heteroclinic orbits generated by the time delay 
term, uncovered by numerically integrating the set of equations (JT|) using a fourth order 
Runge-Kutta method adapted to the time-delay approach. 

As in the case of the FP analysis, a detailed investigation of all possible trajectories 
of the system ([1]) for large regions of parameter space has been carried out. For the sake 
of clarity, here we report only the results that illustrate the typical patterns obtained in 
different regions of the parameter space. Since FP is unstable, almost all trajectories 
starting in its neighborhood are pushed away. When r = the system evolves rapidly 
to the infected state {FPi) and the onset of AIDS occurs a few weeks after the infection. 
This would correspond to cases where the disease progresses rapidly , which in general 
is observed in HIV patients who have been contaminated under immune-depression 
conditions, caused by malnutrition or other diseases that compromise TCDA + cells 
(e.g. tuberculosis). If we increase the value of r from to 4, the trajectory is deviated 
from FP% towards the neighborhood of FP , where the slow transit leads to the a steady 
increase of latency period. 

Figure 2(a) summarizes the main results obtained for this model: the appropriate 
description of the three-stage dynamics observed in infected individuals who have not 
been submitted to drug therapies. In contrast to other ODE models, here, using a 
unique set of parameters, we were able to reproduce the primary infection and a large 
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Figure 2. (a) Time evolution of HIV infection showing three distinct phases. 
Squares, circles and triangles represent h\ + h 2 , A + B and d cells, respectively. The 
parameter values and initial conditions adopted are: /ci=0.163, A;2=0.228, fc3=0.65, 
fc 4 =3.25xl(r 5 , fc 5 =0.606, fc 6 ^0.02, n=4, p=l, g=1.13, /ii=0.95, A^O.Ob. (b) The 
same orbit in a plane projection of phase space. Diamond indicates the initial 
point. Inset: details of the first transit close to FP . The trajectory first lies in 
the attracting manifold, but later it is pushed away from FPq along the repulsive 
manifold. Hcteroclinic orbits occur if FP is reached. 
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range of latency periods (plateau in Figure 2a) that vary from weeks to years. Moreover, 
a plateau is obtained on the T cell counts when the trajectory on the phase space goes 
towards the region of very slow dynamics, close to FPq (see Figure 2b and inset). As 
previously mentioned, the deviation of the trajectory in the phase space to the slow 
transit region largely depends on the values of q and r. Since an increase of these 
parameters leads to a decrease in i?e(A 2 ,3) , the return of the phase space trajectory 
towards FPi proceeds at a much slower pace, thus increasing the length of the latency 
period. 

With regard to the estimates of the latency period, we observe that at the end of 
the plateau, the healthy TCD4 + cell counts oscillates for some time until it stabilizes 
below the threshold associated to the onset of AIDS. Since the latency period is defined 
as the time between the end of the primary infection and the onset of AIDS, we may 
make a second estimate of the latency period that will be larger than that roughly 
estimated from the plateau. The parameter values used in Figure 2 were chosen so as 
to maximize the length of the latency phase. However, the same qualitative pattern, 
with a well-defined latency period, is observed for a finite region of the parameter space. 
Usually, the three-phase pattern remains stable to changes ~ 1 — 2% in parameter values 
of the reported sets. Besides this, if no distinction is made between hi and h 2 , (p = q), 
similar patterns are obtained, albeit with a smaller plateau. For other regions of the 
space parameter, the decay of the trajectory to FP\ may proceed with less defined slow 
transits close to FP and therefore lose the clearly defined latency period. 

Different patterns other than those of the three-stage dynamics may be obtained for 
other parameter values not considered here. Without going into detailed description, a 
brief comment follows on the most important patterns obtained by changing the values 
of q and r. For low values of q, the plateau is reduced, however, if it is increased beyond 
1.18, we observe that all trajectories starting in a neighborhood of FP converge to 
FPq. These heteroclinic trajectories, linking two unstable FPs that are not related to 
the infected state, describe the situation well documented in the literature, where the 
patient has contact with the virus but does not develop the disease. Moreover, for 
q = 1.13 and r = 3, the latency phase is still short, while, for r = 5, a sustained 
oscillation pattern is found due to the existence of stable limit cycles around a locally 
stable FP\. In such a case, two attracting stable sets are present. In other words for 
r > 4.5 the monotonic behavior disappears. 

Therefore, it is possible 'grosso modo' , to group the different trajectories into the 
following groups: i) rapid decay to FP\\ ii) decay to FP\ through a latency phase, whose 
duration will depend on the values of q and r; iii) trajectories converging to FPq for 
q sufficiently large (>1.18); iv) stable limit cycle around FPi with a large amplitude, 
when r is larger than 4.5 

Figure 3 illustrates the average values of the cell counts resulting from the 
integration of ([1]) for N t = 200 trajectories. The initial conditions were randomly 
chosen from the very close neighborhood of FPq for the same parameter values used in 
Figure 2. This assumption may be interpreted as if each of the 200 individuals had been 
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Figure 3. Average values of healthy, infected and dead cell counts resulting from 200 
trajectories using different initial conditions starting in a neighborhood of FPq. The 
parameter values are the same as those used in Figure 2, as are the symbols used to 
describe different concentrations . The error bars indicate the standard deviation. 

exposed to different initial viral loads. Similar results are obtained if slight fluctuations 
of the parameters are permitted to characterize different individuals of the simulated 
population. When we compare the results obtained to its corresponding counterparts 
obtained with the CA model [15] , we observe the same overall dispersion pattern in both. 
However, in the present case the average dispersion is lower in the primary infection and 
after the onset of AIDS, but as obtained with the CA model it increases greatly during 
the latency phase. The large error bars obtained during the latency period indicate a 
large range of variability on the development of the disease among the 200 individuals 
(trajectories). 

The results of this study show that a set of ordinary differential equations with time 
delay terms describing the retard on mounting the specific response, is able to provide 
a description of the entire course of the HIV infection using a single parameter set. To 
our knowledge, this has not been achieved previously by any other similar approach. 
As the model is essentially based on the rules of a CA model, our results indicate that 
these rules capture the main elements and mechanisms underlying the dynamics of the 
infection of TCDA + cells. In the CA model, spatial localization is very important for 
the emergence of distinct time-scales. In the model presented here, with no spatial 
dependence, the individual phases are found to be associated with particular structures 
in the phase space, which are induced by the time-delay terms. These terms describe the 
interplay amongst new generations of infected cells (strains) and new specific responses 
and increase the number of possible classes of solution. Although the number of different 
solutions is not as broad as those of systems with explicit spatial dependence, it is 
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sufficient to generate phase space structures that drive the system to regions of very- 
slow dynamics. A more detailed discussion of all types of solution for system ([1]) will 
be published elsewhere. 
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